
;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   In this example we first get the entire field over time, which will
;   make it easier to calculate tendencies

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
;  a = addfile("../wrfout_d01_2000-01-24_12:00:00.nc","r")
;  filename = systemfunc("ls /public/users/wrf/wrfout/2012060100/wrfout_d01_*")

  if (.not. isvar("input_file")) then
      if (.not. isvar("start_date")) then
          start_date = systemfunc("date +%Y%m%d") + "00"
      end if
      input_file = systemfunc("ls /public/users/wrf/wrfout/" + start_date + "/wrfout_d01_*")
  end if

  if (.not. isvar("output_prefix")) then
      output_prefix = "/public/users/wrf/wrfpost/test/precip24h/precip24h_d01_" + start_date + "_"
  end if

  print(input_file)

  a = addfile(input_file + ".nc", "r")


; We generate plots, but what kind do we prefer?
;  type = "x11"
type = "pdf"
; type = "ps"
; type = "ncgm"

  wks_type = "png"
  wks_type@wkWidth = 800
  wks_type@wkHeight = 800


; Set some basic resources
  res = True
  res@MainTitle = "WRF Precipitation"

  pltres = True
  pltres@NoTitle = True
  pltres@CommonTitle = True
  pltres@NoHeaderFooter       = True

  mpres                             = True
  mpres@NoHeaderFooter              = True
  mpres@mpDataBaseVersion           = "MediumRes"
  mpres@mpDataSetName               = "Earth..4"
  mpres@mpGeophysicalLineColor      = "Black"
  mpres@mpNationalLineColor         = "Black"
  mpres@mpUSStateLineColor          = "Black"
  mpres@mpGridLineColor             = "Black"
  mpres@mpLimbLineColor             = "Black"
  mpres@mpPerimLineColor            = "Black"
  mpres@mpGeophysicalLineThicknessF = 2.0
  mpres@mpGridLineThicknessF        = 1.0
  mpres@mpLimbLineThicknessF        = 1.0
  mpres@mpNationalLineThicknessF    = 2.0
  mpres@mpUSStateLineThicknessF     = 1.0
  mpres@mpOutlineBoundarySets       = (/"National","Continental", "Geophysical"/)   ; turn on country boundaries
  mpres@mpOutlineSpecifiers         = "China:states"
  mpres@mpGeophysicalLineColor      = "Black"       ; color of cont. outlines
  mpres@mpGeophysicalLineThicknessF = 2          ; thickness of outlines
  mpres@mpNationalLineThicknessF    = 2          ; thickness of outlines



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Reading times
  char_times = a->Times
  hours = WRF_Times2Udunits_c(char_times, 0)
  ntimes = dimsizes(hours)         ; number of times in the file
  wrf_dt = (hours(ntimes-1) - hours(0))/(ntimes-1)
  ; convert to Beijing Time
  hours = hours + 8
  times_cst = ut_string(hours, "%Y-%N-%D_%H:%M:%S")
  timesTideStr_cst = ut_string(hours, "%Y%N%D%H%M%S")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  precip24_dt = 24
  dt = tointeger(precip24_dt / wrf_dt)
  do it = dt,ntimes-1,dt             ; TIME LOOP - start at hour 3 as we interested in 3hourly tendencies

    output_file = output_prefix + timesTideStr_cst(it)
    wks = gsn_open_wks(wks_type, output_file)

    print("Working on time: " + times_cst(it) )
    res@TimeLabel = times_cst(it)   ; Set Valid time to use on plots
    mpres@NoHeaderFooter       = True

  ; First get the variables we will need        

    slp = wrf_user_getvar(a,"slp",it)    ; slp
      wrf_smooth_2d( slp, 3 )            ; smooth slp

    tc2 = wrf_user_getvar(a,"T2",it)     ; T2 in Kelvin
       tc2 = tc2-273.16                  ; T2 in C
    tc2@units = "C"
      wrf_smooth_2d( tc2, 3 )            ; smooth slp
    

  ; Get non-convective, convective 
  ; Calculate total precipitation
    rain_exp = wrf_user_getvar(a,"RAINNC",it)
    rain_con = wrf_user_getvar(a,"RAINC",it)
    rain_tot = rain_exp + rain_con
      rain_tot@description = "Total Precipitation"

    rain_exp_tend = wrf_user_getvar(a,"RAINNC",it) - wrf_user_getvar(a, "RAINNC", it-dt)
    rain_con_tend = wrf_user_getvar(a,"RAINC",it) - wrf_user_getvar(a, "RAINC", it-dt)
    rain_tot_tend = rain_exp_tend + rain_con_tend
      rain_exp_tend@description = "Explicit Precipitation Tendency"
      rain_con_tend@description = "Param  Precipitation Tendency"
      rain_tot_tend@description = "Precipitation Tendency"


   ; Plotting options for Sea Level Pressure
     opts_psl = res          
     opts_psl@ContourParameters = (/ 900., 1100., 4. /)
     opts_psl@cnLineColor       = "Blue"
     opts_psl@cnInfoLabelOn     = False
     opts_psl@cnHighLabelsOn    = True
     opts_psl@cnLowLabelsOn     = True
     opts_psl@cnLineLabelPerimOn = True
     opts_psl@cnLineLabelFontHeightF = 0.02
     opts_psl@gsnContourLineThicknessesScale = 1.5
     opts_psl@NoHeaderFooter       = True
     contour_psl = wrf_contour(a,wks,slp,opts_psl)
     delete(opts_psl)
    

   ; Plotting options for Precipitation
     opts_r = res
     opts_r@NoHeaderFooter       = True
     opts_r@UnitLabel            = "mm"
     opts_r@cnLevelSelectionMode = "ExplicitLevels"
     opts_r@cnLevels             = (/.5, 1.0, 5.0, 7.0, 10.0, 20.0, 25.0, 50.0, \
                                     70.0, 100.0, 250.0/)
     colorMap                    = (/ (/255,255,255/), (/197,237,4/), (/146,221,7/), \
                                    (/102,208,9/), (/30,193,14/), (/23,183,18/), \
                                    (/23,170, 20/), (/50,160,80/), (/255,255,0/), \
                                    (/255,165,0/), (/255,0,0/), (/208,32,144/) /) / 255.0
     opts_r@cnFillColors         = colorMap
     opts_r@cnInfoLabelOn        = False
     opts_r@cnConstFLabelOn      = False
     opts_r@cnFillOn             = True

   ; Total Precipitation (color fill)
     contour_tot = wrf_contour(a,wks, rain_tot, opts_r)
 
   ; Precipitation Tendencies 
     opts_r@SubFieldTitle = "from " + times_cst(it-dt) + " to " + times_cst(it)
 
     contour_tend = wrf_contour(a,wks, rain_tot_tend,opts_r) ; total (color)
     contour_res = wrf_contour(a,wks,rain_exp_tend,opts_r)   ; exp (color)
     opts_r@cnFillOn = False
     opts_r@cnLineColor = "Red4"
     contour_prm = wrf_contour(a,wks,rain_con_tend,opts_r)   ; con (red lines)
     delete(opts_r)

    ; Plotting options for T
      opts = res                         
      opts@NoHeaderFooter       = True
      opts@cnFillOn = False
      opts@cnLineColor = "Red4"
      contour_tc = wrf_contour(a,wks,tc2,opts)
      delete(opts)


   ; MAKE PLOTS

     ; Total Precipitation Tendency + SLP
       plot = wrf_map_overlays(a,wks,(/contour_tend, contour_psl/),pltres,mpres)

    delete(wks)
  end do        ; END OF TIME LOOP

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

end
